###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*##
###*### COMPUTE LOCAL COUNTERFACTUAL EXPERIMENTS          #*###*###*###*###*###*###*###*##
###*### USING PREVIOUSLY ESTIMATED PARAMETERS             #*###*###*###*###*###*###*###*##
###*### LOCATION 1                                        #*###*###*###*###*###*###*###*##
###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*##

#' FROM THE PAPER:
#' Child Care Markets, Parental Labor Supply, and Child Development
#' BY:
#' Samuel Berlinski, Maria Marta Ferreyra, Luca Flabbi and Juan David Martin (2022)

# ___________________________________________________________________________________ ####
# START ####

library(lubridate)
library(snowfall)
library(readxl)
library(plyr)

rm(list = ls())
gc()

## Directory
folderpath <- '~/Documents/BM-BID/Codes/Code_LocalMarkets/'
setwd(folderpath)

# ___________________________________________________________________________________ ####
# INPUTS ####

## Functions
source('functions.R')
source('simul_market.R')
source('simul_market_PP.R')

## Sim data
load('SimData_L1.RData')

## Estimates
load('EstimatesL1_01.RData')

# ___________________________________________________________________________________ ####
# SETUP ####

## Market given settings
J     <- nrow(data)
NL    <- 5
ratio <- 0.9062
minw  <- 5.15*ratio
M     <- 40*NL/J/(0.3 * 1.25); M
Fixed <- c(128.90, 322.25)*ratio
Cn    <- minw
Cr    <- 0
I     <- minw*8
subs  <- 0.35
Nrel  <- round(J*0.4*1.0776805)
mini  <- 80

## Observed heterogeneity
q0  <- data$q0
ed  <- data$edu
edf <- data$edf
wf  <- data$wf
sm  <- data$single
To  <- 16 + 2*sm

## Dad's wage implied random term
ef  <- log(wf)      
mwf <- mean(ef[sm == 0])
swf <- sd(ef[sm == 0])
ef  <- (ef - mwf)/swf
ef[wf == 0] <- 0

## Random draws
set.seed(123)
ev   <- matrix(rnorm(3*J), J, 3)
ew   <- matrix(rnorm(J),   J, 1)
ec   <- rowMeans(apply(matrix(rnorm(NL*10000), NL, 10000), 2, sort))
ec   <- (ec - mean(ec))/sd(ec)
reli <- sample(1:J, size = Nrel, replace = F)

## Standarization for wages draws
ew[ed == 1] <- (ew[ed == 1] - mean(ew[ed == 1]))/sd(ew[ed == 1])
ew[ed == 2] <- (ew[ed == 2] - mean(ew[ed == 2]))/sd(ew[ed == 2])
ew[ed == 3] <- (ew[ed == 3] - mean(ew[ed == 3]))/sd(ew[ed == 3])

# ___________________________________________________________________________________ ####
# SIMULATE BL ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1)
sfStop()

## Save
market_BL <- market
save(market_BL, file = 'simL1_BL.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF1: CASH TRANSFER ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1, Cash = 40)
sfStop()

## Save
market_CF1 <- market
save(market_CF1, file = 'simL1_CF1.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF2: VOUCHER UNCOND. ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1, tau = 5)
sfStop()

## Save
market_CF2 <- market
save(market_CF2, file = 'simL1_CF2.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF3: VOUCHER MOM WORK ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1, tau = 5, both = T)
sfStop()

## Save
market_CF3 <- market
save(market_CF3, file = 'simL1_CF3.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF4: VOUCHER INCOME ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1, tau = 5, both = F, mini = mini)
sfStop()

## Save
market_CF4 <- market
save(market_CF4, file = 'simL1_CF4.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF5: VOUCHER QUALITY ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1, tau = 5, hqv = T)
sfStop()

## Save
market_CF5 <- market
save(market_CF5, file = 'simL1_CF5.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF6: PUBLIC PROVISION P = MeC ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market_PP(pp1)
sfStop()

## Save
market_CF6 <- market
save(market_CF6, file = 'simL1_CF6.RData')

# ___________________________________________________________________________________ ####
# END ####